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Abstract 

The critical behavior of the Ising model with non-conserved dynamics and an external driv- 
ing field mimicking a shear profile is analyzed by studying its dynamical evolution in the short 
time regime. Starting from high temperature disordered configurations (FDC), the critical 
temperature T c is determined when the order parameter, defined as the absolute value of the 
transversal spin profile, exhibits a power-law behavior with an exponent that is a combination 
of some of the critical exponents of the transition. For each value of the shear field magnitude, 
labeled as 7, T c has been estimated and two stages have been found: 1) a growing stage at 
low values of 7, where T c ~ 7* and ip = 0.52(3); 2) a saturation regime at large 7. The same 
values of T c {pj) were found studying the dynamical evolution from the ground state configura- 
tion (GSC) with all spins pointing in the same direction. By combining the exponents of the 
corresponding power laws obtained from each initial configuration the set of critical exponents 
was calculated. These values, at large external field magnitude, define a new critical behavior 
different from that of the Ising model and of other driven lattice gases. 

PACS numbers: 75.30.Kz; 64.75. +g; 05.45.pq; 05.70.Ln. 

1 Introduction 

The statistical mechanics of equilibrium phenomena is a very useful theoretical framework for un- 
derstanding the thermodynamic properties of many-particle systems from a microscopical point 
of view. However, in nature, most of the systems evolve under out-of-equilibrium conditions, and 
there is not yet a suitable general framework to study them as in the case of equilibrium systems. 
Nevertheless, some progress have been achieved in the knowledge of far from equilibrium behavior 



by means of simple models, capable to catch the essential physics of non-equilibrium processes. 
In this context, we introduce a very simple model, derived from the Ising model, driven out of 
equilibrium by an external field that mimics the effects of a uniform shear profile pQ . This model 
evolves with a non-conserved dynamics, corresponding to model A in the classification of Hohenberg 
and Halperin [2], and it was already used by Cavagna et. al [3] and by Cirillo et. al [4j to study 
phase separation. 

We will focus on the study of phase transition properties in this model. Typical configurations 
observed in our simulations are displayed in figure [TJ At low temperatures the system appears 
ordered with elongated domains directed along the field direction. At high temperatures the system 
exhibits a gas-like appearance with disordered patterns. Similar ordered and disordered phases, also 
experimentally found [5], generally characterize the behavior of sheared binary systems. As usually 
in systems with an applied external field, the transition point is a function of the magnitude of the 
driving field [6]. 




Figure 1: Snapshot configurations corresponding to the two phases of the Ising model with shear. 
On the left panel we observe a typical stripe-like configuration at T < T c . On the right panel a 
configuration at T > T c is displayed. The external field is applied along the horizontal axis. 

Previous theoretical and experimental studies have shown that sheared binary systems undergo 
a second order phase transition at a critical temperature T c [6]. In diffusive systems the effect of 
the external driving field is to inhibit fluctuations so that the critical temperature is expected to 
increase with the magnitude of the driving. In a continuum model with non-conserved dynamics, in 
the large- N analytical approximation, it has been found that the value of the critical temperature 
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depends on the driving field following a power law at small field magnitudes [7] . In previous Monte 
Carlo studies on the critical behavior of sheared Ising models [8], it was not possible to extract 
information about the critical temperature, due to numerical uncertainties and finite size effects [9]. 
In view of this, we revisit this issue, in order to determine the critical temperature of the model 
as a function of the magnitude of the external field, and to compute for the first time the critical 
exponents in this model. For the sake of comparison, we will contrast the obtained values with 
those computed for the 2d driven lattice gas model (DLG) pTOj that will be briefly described in the 
next section. 

To study the phase transition in the model, the critical dynamical behavior will be investigated 
by monitoring the time evolution of some observables before the system reaches nonequilibrium 
steady states (NESS). This technique, generally called short time dynamics fT2] is an alter- 
native convenient way to obtain both the critical temperature and exponents precisely, with less 
computational cost than other methods commonly used, such as the finite size scaling applied to 
the specific heat and response functions, that require a considerably amount of simulation time 
in order to reach NESS. Furthermore, since the measurements are carried out in the first steps of 
evolution, the short time dynamic approach is free of the critical slowing down. 

The manuscript is organized as follows: In Section [2j the Ising model with the external shear 
field is introduced. The technique used to study the phase transition is described in Section 03 The 
simulation results are presented in Section HI and finally the conclusions are stated in Section \E[ 

2 The model 

We will consider the nearest-neighbor two-dimensional Ising model with a single-spin-flip ther- 
malization dynamics, e.g. the Metropolis dynamics [13]. The driving field will be defined in order 
to mimic the convective velocity shear profile 

vx(y) = iy v y = (l) 

where the parameter 7 is called the shear rate and represents the shear field magnitude. If the 
system is imagined as a sequence of layers labelled by y, then jy is the displacement of the layer 
y in a unit of time. If L y is the vertical size and %ax is the speed of the fastest layer, then 
jL y = Umax- 

The model is defined on a square lattice A of horizontal and vertical sizes L x , L y respectively, 
with periodic boundary conditions in the L x direction and free in the L y direction. More precisely, 
let f2 = { — 1, +1} A be the space of configurations and, for a S fi, let a x ^ y be the value of the spin 
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associated to the site (x, y) € A. Then the Hamiltonian of the model is 

(2) 

y— lx- 1 x — 1 y—1 

with aL x +i, y = <Ti, y for all y = 1, . . . , and J is a positive real coupling constant, which means 
that the interactions are ferromagnetic. We will combine the thermalization dynamics with an 
algorithm introducing the shear in the system. The shear is superimposed to the thermalization 
dynamics with typical rates not depending on the thermalization phenomenon, but fixed a priori. 
This has been done in different ways in [31 [El [14]. In this paper we use a very ductile generalization 
of those dynamics aiming to introduce the shear effects in a way resulting competitive with respect 
to the thermalization process. Notice that our dynamics results from the combination of two steps: 
i) a thermalization step which would bring the system in the usual equilibrium; ii) a shear step 
which changes the configurations of the system forbidding to reach the equilibrium. Therefore, 
all together, our algorithm does not satisfy local detailed balance expressed in terms of standard 
equilibrium probabilities of configurations. Similar models have been also used in different context 
of non-equilibrium studies [15J. 

Let the time unit be the time needed for a full thermal update of the entire lattice, e.g. a full 
sweep of the Metropolis algorithm. The shear algorithm is parametrized with a submultiple r of 
L y L x (the period of the shear procedure), a positive integer A < L x /2 (the number of unit cells that 
a row is shifted when the shear is performed), and a non-negative real v < 1/L y . The dynamics of 
the model that we study in this paper is defined in a precise way via the following algorithm: 

1. set t = 0, choose <jq e il, and set n = 0; 

2. increase by 1 the index n, and choose at random with uniform probability l/L x L y a site of 
the lattice and perform the elementary single-site step of the thermalization dynamics; 

3. if n is multiple of r a layer is randomly chosen with uniform probability l/L y . Then, if y is 
the chosen layer, all the layers with y > y are shifted by A lattice spacings to the right with 
probability vL y \ 

4. if n < L x L y goto [21 else denote by at+i the configuration of the system; 

5. set t = t + 1, set n = 0, and gotoO 

We note that if v — X/L y the shift at step [3j is surely performed and this case will be later addressed 
to as full shear. The smoothness of the shear field, eqs. (0}, is ensured by the random choice of the 
layer y in step 3. Now we want to express the shear rate 7, introduced in equation ([l]), in terms 
of the parameters of our dynamics. We have to estimate the typical displacement per unit of time 
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of the row labelled by y. Such a row is involved in a shear event, step [3] of the algorithm above, if 
and only if the extracted row y is such as y < y, and this happens with probability y/L y . Since the 
shear event results in a shift with probability vL y , the probability that during a shear event the 
row y does shift is given by 

V 

-j— x vL y = vy. 

Ly 

By noting that the number of shift events per unit of time is equal to L x L v /t and recalling that 
the shift amplitude is A, we have that the typical shift of the row y per unit of time is given by 



T 

By using definition (U) we finally get 7 = L x L v v\/t, which becomes 7 = L x X/t in the case of full 
shear. 

It is important to remark that there exists a large variety of models that evolve under non equilib- 
rium states by the action of an external field. An example is the driven lattice gas (DLG) where 
the driving field is not superimposed to the thermalization dynamics, but it is rather inserted in the 
Metropolis transition rates, that become anisotropic, and biases the movement of particles along its 
direction [10J. Furthermore, it exhibits a second order phase transition at particle density p — 1/2, 
between an ordered phase at low temperatures characterized by regions of low and high particle 
density, called stripes, oriented along the field direction, and a disordered phase at high temper- 
atures with the appearance of a lattice gas. Both ordered and disordered phases have a similar 
appearance with those exhibited in figure [U The critical temperature increases with the magnitude 
of the external field, saturating at T c ~ 1. 41X^(0) in the case of infinite [10]. Here, T c (0) — 2.269 
J/ks is the critical temperature of the 2d Ising model (J is the coupling constant and ks is the 
Boltzmann constant). 



3 Dynamical Critical Behavior 

It is known that, for systems exhibiting critical behavior, the relevant observables measured in 
equilibrium stationary states can be written in terms of power laws, with characteristic critical 
exponents due to the divergence of both the spatial correlations and the correlation time. In recent 
years, however, the attention has been also focused to the early stages of the evolution of the sys- 
tem towards the critical state, that is, to a microscopic time regime where the spatial correlation 
length is small compared with the system size [TT]. Within this regime, it is possible to measure 
scaling laws of the observable quantities OUS]. This new method to study second-order phase 
transitions, called short time dynamics, allows to estimate the critical temperature and to compute 
the critical exponents of the transition with relative quickness and avoids the shortcomings that 
more usual techniques present to study critical behavior. Furthermore, the short time dynamics 
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has been applied to investigate the critical behavior of a wide range of systems of different nature, 
such as models showing criticality under equilibrium conditions, such as e.g, the XY systems [17], 
the 2d 3-state Potts model [l8| . the Ising magnet under different lattice geometries OEOj, and 
of nonequilibrium critical models such as the driven diffusive lattice gas (DLG) [IH E21 EH] , etc. 
In these last three works, a detailed analysis of the second-order |21l [22] and first-order [23] non- 
equilibrium phase transition was performed by using the short time critical dynamic methodology. 
For the second-order phase transition, the excellent agreement between critical exponents evaluated 
using the standard (stationary) and dynamical (short time) approaches strongly support the ro- 
bustness of this method. Encouraged by this success, our goal is to extend the short time dynamics 
concept to this model, basing our ideas on the already developed short time dynamics method for 
the DLG model in ref. [21]. 

The above mentioned scaling laws can be observed employing two different initial configurations, 
namely: 1) Fully disordered configurations (FDC's), which means that the system is initially placed 
in a thermal bath at T — > oo , and the system configuration is similar to that exhibited on the right 
panel of figure [TJ 2) Completely ordered configurations or ground state configurations (GSC's) as 
expected for T = 0. In our model, based on the fact that the equilibrium Ising model has all the 
spins pointing in the same direction (i. e. magnetization M = 1 or -1) at this temperature, we will 
adopt this configuration as the ground state for testing the short time dynamic behavior. 
The shear field introduces anisotropic effects, that generates anisotropic correlations in the system. 
As a consequence of this, there will be two correlation lengths, namely: 1) A parallel or longitudinal 
correlation length £|| along the external field direction, and 2) A perpendicular or transverse cor- 
relation length £_l perpendicular to £y. Whatever the initial condition is used to start the system, 
both spatial correlation lengths are quite small or zero at the beginning of the dynamic process, and 
near the critical temperature T c they increase dynamically as a power law £||(±) x t 1 ' z ^ ± '> , where 
2|i(_L) is the dynamic critical exponent in the respective directions. 

Before we start to describe the theoretical basis of the technique applied to this model, we set the 
external driving field along the horizontal direction, i.e. the x axis. Also, we need to define quanti- 
ties that are relevant in the critical behavior of the model. Based on the morphological appearance 
of typical configurations present in the system (see figure [I]), we will consider a variant of the order 
parameter OP employed in the critical study of the DLG model [21] : 



The order parameter defined in this way can take into account the small ordering that appears at 
the early stages of the evolution. 

There is one more point to take into account before we start to expose the method applied to 




(3) 



where P(y) — -jj- 53a™i a x,y ls the average of the spin profile in the shear field axis. 
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this system. In all formulas below we will assume, and demonstrate later, that only the parallel 
correlation length is relevant in the short time critical evolution of the system. In fact at T ~ T^j) 
and at early times of evolution, parallel and perpendicular correlations begin to increase. However, 
domains of perpendicularly correlated spins are broken by the shear, and assume a characteristic 
elongated shape, also observed in many experimental studies of sheared systems. As a consequence 
of this, transversal correlations grow slower than parallel correlations, so they do not take part in 
the dynamic critical behavior of the model at short times. This effect was also shown for the DLG 
model (211 [22] . Since this happens independently of the initial configuration, we will take z = z\\ in 
every expression below. Furthermore, they must contain the anisotropic finite size dependence in 
order to match the usual anisotropic scaling forms for the NESS regime. 
Starting with FDC's, the scaling law proposed for the order parameter (OP) reads |21j . 

OP(t,<l>,L x ,L y ) = b-^»WP*(b- z t,b 1 ^cb,b- 1 L x ,b~^^L y ), (4) 

where t is the time, b is the spatial rescaling factor, (3 is the critical exponent of the order parameter, 
are the correlation length critical exponents in the x(y) axis oc (/) _l/ ii( J -> ), OP* is a scaling 

function, z is the already mentioned dynamic exponent of the longitudinal correlation length, and 
= T ~J C . Notice also that L y is rescaled by b^ u± ^« to include possible shape effects in the 
dynamic critical behavior |10j . 

To generate the FDC initial conditions, the lattice is filled at random with exactly poL x L y particles, 
being po = 1/2 the density of up spins. However, the number of particles on each row parallel to the 
field axis is not the same for all rows. This generates tiny density fluctuations along this direction, 
which are of the order of (1/L X )~?, in agreement with the central limit theorem. According to 
equation ([4]), these fluctuations add up, and the amplitude of OP depends on L x We have to 
take into account this expression for the final form of the time evolution of OP. Setting b ~ t x / z , 
eq. |(4j) becomes: 

OP(t,^,L x ,L y ) = t-Pl v \\*OP\\,t x l v \\ z ^t- xlz L x ,t- v ±l v \\ z L y ). (5) 
Then, if t~ x l z L x is extracted out of the scale function in Eq. ©, we have: 

OP{t, <f>, L x ,L y ) = t- f3 ^ z (t-^ z L x ) x OP**(l, t 1 ^, t- u ^ z L y ), (6) 

1 

but since OP ~ 1/LJ , then x — —1/2, so, the final expression for OP(t) is the following: 

OP(t, 4>, L x ) = L-VH^OP**^*® L x , L y oo, (7) 

1 These fluctuations are equivalent to the initial magnetization mo in the original formulation of short time 
dynamics [16] . 
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with ci = (1 - 2(3/v\\)/2z [2T] . 



Furthermore, it is easy to show that the logarithmic derivative of OP with respect to 0, given 
by Eq. J7]) at criticality, behaves as 

dj,lnOP(t,cj>) (xt c \ (8) 

where the exponent is ci — \jv\\z. 

On the other hand, starting the system from the GSC configuration described above, and according 
to the scaling behavior proposed in [21J, we have 

OP(t, 0, L x , L y ) = b-^^OP***(b- z t, b- x L x ,b- v wl v ±Ly), (9) 

where OP*** is another scaling function . Here we have also included the shape scaling factor 
lr v - s -l v \\L y . Proceeding in the same way as in the above case, we have, taking b ~ t x / z in (|9|) at 
criticality, the following expression for OP: 

OP(t)oct~ C3 L x ,L y ->oo, (10) 

with an exponent C3 = f3/v±z. Moreover, the derivative of Eq. IjlOp with respect to 4> at criticality 
is given by 

d^OP(t) o,t Ci , (11) 
where the exponent is C4 = (1 — j3)/v±z. 



4 Simulation Results 
4.1 Critical Temperature 

In this work we used rectangular and square lattices of different sizes L x ,L y , in the range 128 < 
L X i L y < 10000 lattice units. The critical dynamics of the model was investigated as a function of 
the shear magnitude, in the range 1/32 < 7 < 50. The temperatures were measured in units of 
J/ks, being the Boltzmann constant, and the time is measured in Monte Carlo steps (mcs), 
where one unit consists of L x L y attempts for spin updates. The time evolution was sampled from 
100 to 1000 realizations of the system, according to each initial condition and temperature. 
We begin showing our results by considering the critical dynamic evolution of the system when it is 
started from FDC configurations, and then coupled with a thermal bath at T ~ T c . In figure [21 the 
time evolution of OP(t) is shown for a system where a shear field is applied with shear rate 7 = 5. 
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The best power law behavior is obtained at T = T c = 2.660 before OP(t) reaches a saturation value 
due to finite size effects. On the other hand, if T < T c or T > T c , OP deviates from the power law 
behavior as equation (0 states, and the curve shows an upward or downward bending, respectively 
(see the curves corresponding to T — 2.655 and T = 2.665 in the same figure) 




Figure 2: Log-log plot of the short time evolution of OP(t) starting from FDC configurations, in a 
system with L x = L y = 512 and 7 = 5. The best fit of the raw data gives a power law behavior at 
T c = 2.660, as it is indicated by the dashed straight line. Upward and downwards deviations from 
this behavior can also be observed for T — 2.655 and T = 2.665, respectively. 

The critical dynamic evolution was investigated by performing simulations also on rectangular 
lattices. The plots in figures [3] display the dynamic evolutions of OP(t) at T = T c for two shear 
field magnitudes, 7 = 1/2 and 7 = 10 in a) and b) panels, respectively. The lattice sizes used are 
indicated in the legend of each plot by the notation L x x L y . In the main plots of each figure, all 
early-time evolution exhibits the same power-law behavior with similar values of the exponent c± 
(equation Q). Then, a saturation value is reached, OP sa t, that depends only on L x , as it can be 
observed in lattices with longitudinal sizes L x — 500 and L x — 1000 for 7 = 1/2 and L x = 500 for 
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7 = 10 respectively. So, these plots show that the early-time critical evolution of the system is free 
from lattice shape effects [lOj, because the relation between the longitudinal and transversal sizes 
is different for each studied case. The same was observed in the short-time critical evolution of the 

1 /2 

DLG model [21]. In addition, the power law behaviors can be collapsed by rescaling OP(t) by Lj 
as it is proposed in equation J7]) . This is shown in the insets of both figures 
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Figure 3: Log-log plots of the time evolutions of OP{t) at T = T c , for the system in rectangular and 
square lattices of sizes L x x L y , as indicated in the legends. In a) the external field has magnitude 
7=1/2 while in b) the shear field magnitude is 7 = 10. The straight lines are least-square fits of 

the data. The insets in each figure show the collapse of the power law behaviors when OP(t) is 

1/2 

multiplied by L x 

Then, the critical points of the system at several values of the shear field magnitudes 7 were 
found. Figure H] shows that the power-law behavior of OP collapses for large values of 7, i.e. 
7 = 5, 10 and 50, and occurs at approximately the same critical temperatures for each shear rate. 
These temperatures are larger than the estimated for the equilibrium Ising model, T c (j = 0) = 2.269 
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(see table 1). However, the time evolution of OP depends on the shear field if 7 is small, as it can 
be seen for the cases 7 = 1/32 and 7 = 1/10. Although these magnitudes are quite small, they are 
enough to raise the critical temperature to T c (j — 1/32) = 2.29 and T c (j — 1/10) = 2.395, both 
close but greater than the critical temperature of the Ising model. This confirms that the critical 
temperature depends on the shear field magnitude, as found by theoretical studies [7]. 




Figure 4: Log-log plots of the time evolutions of OP(t) at T = T c corresponding to systems with 
different shear fields. The magnitudes 7' are indicated in the legend. 

Once that the critical temperatures T c corresponding to the different 7's were collected, a di- 
agram of critical temperatures versus 7 can be performed. Figure [5] shows that two regimes can 
be distinguished. In the first regime, the critical temperature grows with 7 as a power law, i. e. 
T c (-j)/T c (0) — 1 oc 7*\ The value of the exponent was estimated in V = 0.52(3), which is con- 
sistent with that calculated theoretically in [7j- I n this work (ref. [7]), the critical transition in 
the if 3 — ► (if 2 ) f approximation was studied in a scalar field model based on a convection-diffusion 
equation with Landau-Ginzburg free energy, with the average (if 2 ) self-consistently determined. 
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Above the lower critical dimension, the exponent ip was evaluated to be 1/2 and 1/4 for the cases 
with non-conserved and conserved order parameter respectively. 

Then, T c crosses over to a saturation regime at larger 7's, saturating at T c (j) ^ 1.18T C (0), 
where T c (0) !~ 2.269 is the critical temperature of the 2D Ising model. The increase of the critical 
temperature by action of the external field, and a posterior saturation regime was also observed in 
the DLG model [10]. We want to remark also that in real fluids there is also a negative contribution 
to the shift of the critical temperature coming from hydrodynamics interactions [24] ■ The total shift 
of T c results to be negative for fluids with low molecular weight [24] , was also found in experiments 
|25| . A review for other systems is given in [6]. 



0.150 



0.100 



0.075 



o 
H 

•|3 0.050 
o 
H 



0.025 




Figure 5: Diagram of critical temperatures T c (j)/T c (0) — 1 versus 7, in log-log scale. The fit of the 
data points in the growing regime gives a slope ip = 0.52(3). 

On the other hand, if the system is started from the GSC initial condition (i.e. magnetization 
equal to 1), and then it is left to evolve at the working temperature T ~ T c , OP{t) decreases, and 
follows a power law behavior at T — T c . Upwards or downwards deviations are observed according 
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if T < T c or T > T c , respectively. Figure [6] exhibits the evolution of OP(t) for the same parameters 
of figure [2l A clear power law behavior can be observed at T — 2.66, which is exactly the same 
temperature found when the system was started from FDC configurations. This is precisely the 
signature of a second-order phase transition in the model. In addition, the expected deviations from 
the power law behavior at T = 2.65 and T = 2.67 are also observed. 




Figure 6: Log-log plot of the time evolutions of OP{t) corresponding to a system set in a square 
lattice of size L x = L y = 512 with shear rate 7 = 5. The working temperatures of the thermal bath 
are indicated in the legend. A power law behavior is observed for T c = 2.660 J/fcs- The dashed 
line is a fit of the numerical data. 

Proceeding in the same way we did for the evolutions started from FDC's, the critical behavior 
of the system was investigated when it is started from GSC configurations in rectangular lattices. 
FiguresOa) and b) show the critical time evolution of OP when the system is initiated from a GSC 
configuration, corresponding to 7 = 1/2 and 7 = 10, respectively. It is important to remark that the 
best power law behavior was obtained when the systems evolved at the same critical temperatures 
found when they were initiated from FDC configurations. According to the results, the transversal 
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size L y does not play a relevant role in the critical behavior of the system as it is shown in the 
evolutions in lattices with L x = 500 (figure a) and b)) and also with L x = 2000 (figure [7] b)), but 
rather the longitudinal size L x is relevant. This is in agreement with the results previously exhibited 
in figure O and allow to conclude that that the critical evolution of the system is independent of 
lattice shape effects for both initial configurations. 





Figure 7: Log-log plots of the time evolutions of OP(t) at T = T c , for the system set in rectangular 
and square lattices of sizes L x x L y , as indicated in the respective legends. In a) the external field 
has magnitude 7 = 1/2 while in b) the shear field magnitude is 7 = 10. The dashed lines displaced 
are least-square fits of the obtained data. 
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4.2 Critical Exponents 

Focusing our attention on the power law behavior at the T — T c , the exponents C\ and c 2 corre- 
sponding to the dynamic evolution of OP (equation |(7j)) and its logarithmic derivative d T lnOP(t, r) 
(equation JSJ) can be estimated. Table Q] enlists the obtained values of T c and the mentioned expo- 
nents for small and large values of 7. The values of c\ and C2, for the case of 7 = 1/32 are different 
compared with the those at larger 7, as it was already observed in figure HI Furthermore, the value 
of ci for 7 = 50 is also slightly different from the corresponding cases of 7 = 5 and 7 = 10. This 
seems to be a rare behavior, since there exists a saturation regime for the critical temperature at 
these values of 7's (figure [5]), and we are induced to think that the dynamic behavior of the system 
is independent of the field magnitude in this limit. Opposite to that, the estimated values of C2 are 
quite similar in the limit of large 7's. 



7 


T c 


Cl 


C2 


1/32 


2.29 


0.180(8) 


0.99(2) 


5 


2.66 


0.239(1) 


0.84(1) 


10 


2.675 


0.238(1) 


0.87(1) 


50 


2.675 


0.224(1) 


0.85(1) 



Table 1: Exponents ci and C2, obtained from FDC initial conditions, corresponding to the values 
of 7 enlisted in the first column. 

Table [2] enlists the values of the exponents C3 and C4 that were obtained from a least-square fits of 
the critical evolution of the system when it is initiated from the GSC configurations (see equations 
fl0|) and (HU). Here, the same situation that happened for ci is found for C3. In fact, the value of 
C3 for 7 = 1/32 is different form the rest of the corresponding values at larger 7's, and the value of 
C3 for 7 = 50 is also different from the values estimated for 7 = 5 and 7 = 10. On the other hand, 
the values of C4 are similar for all the reported cases of 7's. 



7 


T c 


C3 


C4 


1/32 


2.29 


0.076(1) 


0.65(2) 


5 


2.66 


0.401(1) 


0.63(1) 


10 


2.675 


0.405(5) 


0.62(1) 


50 


2.675 


0.360(7) 


0.62(5) 



Table 2: Exponents C3 and C4 obtained from the GSC configuration, corresponding to the values of 
7 enlisted in the first column. 

According to the equations developed in section [31 the critical exponents of the second-order phase 
transition, are obtained by combining the estimated exponents ci, C2, C3, and C4 enlisted above, 
corresponding to each case of 7 investigated. Table [3j resumes the obtained results and includes the 
critical exponents of both the Ising and DLG model for the sake of comparison. Since this issue for 
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the case of the DLG model still remains an open problem, the obtained theoretical values from all 
proposed theories exposed in refs. [10] and [27] are included. 



7 


a 


z 


v\_ 




1/32 


0.105(1) 


1.77(7) 


0.78(1) 


0.57(5) 


5 


0.39(1) 


0.88(1) 


1.10(1) 


1.35(2) 


10 


0.39(1) 


0.86(1) 


1.14(1) 


1.34(1) 


50 


0.37(1) 


0.94(1) 


1.10(1) 


1.30(6) 


Ising 2d 


0.125 


2.16 


1 


1 


DLG(E=50) (ref. [10J) 


1/2 


w4/3 


«3/2 


«l/2 


DLG (E=50) (ref. [27J) 


«0.33 


wl.998 


wl.22 


w0.63 



Table 3: Table of the calculated critical exponents for each case of 7 employed. The critical 
exponents corresponding to both theories of the critical phase transition of the Ising and the DLG 
models are also included for comparison. Since there is no anisotropy in the Ising model v±_ =v«=v 
must be read. The results of this table do not fit with the relation tp = 1/uz suggested in [24]. See, 
however, in Sect. 5, the discussion about our results our small shear rates. 

An overview of this table deserves some comments. First, the calculated value of (3 at 7 = 1/32 is 
close to f3 = 1/8 calculated for the 2d Ising model. This similarity could drive us to think that the 
effects of the shear field are negligible, but this is not the case, as it is evidenced by the values of 
the rest of the critical exponents. In fact, anisotropy effects are important, even if a small external 
field, in this case 7 = 1/32, is applied. The dynamic exponent 2(7 = 1/32) = 1.77(7) indicates 
that the correlation length (in the longitudinal direction, see 14.3ft grows faster with time than the 
corresponding one in the Ising model z = 2.16, and the difference between v±(j — 1/32) = 0.78(1) 
and I'm (7 = 1/32) = 0.57(5) reveals an anisotropic critical behavior even at small shear rate values. 
The situation is different for the critical exponents at the largest values of 7 investigated. The 
values of the exponents are similar between each other, suggesting that the critical behavior does 
not depend of the applied field. This fact is also present in figured! where the critical temperature is 
approximately the same for the largest values of 7 used. Furthermore, we also noticed that u± > u« 
for 7 = 1/32, while it happens the opposite at large 7. At the moment, a reasonable explanation 
for this issue is not possible due to the lack of a theoretical framework about the critical behavior 
of this model. 

To end this section, one final comment is appropriate. In view of the values exposed in table 3, the 
computed critical exponents do not belong to the universality classes of the Ising or of the DLG 
models respectively. This fact is not surprising since, as we have seen, the shear rate affects the 
critical behavior of the model by inducing anisotropic effects in the equilibrium model that changes 
its behavior. The case of the DLG model is different. Although both models have a similar phase 
behavior, the particle dynamics defined for the DLG model conserves the number of particles while 
our model does not. This difference will probably affect the values of the critical exponents, and 
in consequence there is no reason to expect that both models will belong to the same universality 
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class. 



4.3 Longitudinal Correlation 

In section [31 it is assumed that the dynamic increase of the longitudinal correlation length £y and 
the breakage of the corresponding transversal one £j_ at T — T c are due to the anisotropy effects 
induced by the external shear field. This will cause that the short-time critical dynamic evolution 
of the system initiated from either the FDC or GSC configurations will depend only on the dynamic 
critical exponent z\\ = z that describes the critical dynamic increase of 

To show our hypothesis, we performed a scaling of the whole curve of the OP(t) for the cases of 
7=1/2 (small driving field) and 7 = 10 (large driving field). We propose a phenomenological 
scaling in the spirit of the scaling form used by Family and Vicsek to describe the roughness growth 
of interfaces [26] (obviously in a different context not related to ours) . This is given by the following 
expression 

OP = L-^f{t/L%), (12) 

where L x is the longitudinal size, according to the results exhibited in figures [3j and [7] a) and 
b), respectively. The exponent Ui, i =FDC or GSC, is the exponent that has into account the 
finite-size critical behavior of OP, and z is the relevant dynamic critical exponent. The idea of 
this scaling form is simple: if all the curves can be collapsed by using the same dynamic exponent 
z independently of the initial condition used to start the simulations, it can be demonstrated 
numerically that only one correlation is relevant in the critical dynamic behavior. Furthermore, 
equation (fT2|) must contain both the critical dynamic behavior of OP according to equations ((!]) 
and (fTPl at the early times of evolution, and the finite size behavior in the limit of large times 
(t — > 00) where the correlation length is comparable to L x . Therefore we have that f{u) must be 
/ oc (tj ' L z x )~P/ v w z or / cx (t/L z ,)~ l3 / u±z at early times u <C 1, depending if the initial condition is 
FDC or GSC respectively. This fixes the finite-size exponent uii in uofdc — or ^gsc — 
according to the initial configuration used to start the simulations. 

Figures [8] and [9] a) and b) exhibit both the finite-size dependence of OP(t) with the longitudinal size 
L x (insets), and the scaling function f(t/L x ) (main plots), for the small and large external fields, 
represented by 7 = 1/2 (figures [8]) and 7 = 10 (figures [9} respectively. In all plots, the finite-size 
dependence is obtained by calculating the saturated value of OP(t), OP sa t, from figures [3] and [3 
which were plotted versus L x . 
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Figure 8: Log-log plots of the scaling functions f(t/L z ) (see equation Ijl2p ) of the time evolution of 
OP(t) when the system is started from FDC and GSC configurations (main plot of figure a) and b) 
respectively) corresponding to a system with external field of magnitude 7 = 1/2. The insets of both 
figures show the finite-size dependence and the corresponding power-law fit in a double logarithmic 
scale. The values of the corresponding fits are ujfdc — 0.50(1) (inset of a) ) and ujgsc = 0.53(1) 
(inset of b)), respectively. 

As it can be observed in the insets of the plots in figures [8] a) and b) the size dependence of OP 
in the long time regime can be well fitted by a power law as it is proposed in equation ll"12|) . The 
estimated exponents oj/dc = 0.50(1) and Lo gS c — 0.53(1) were not consistent with the expected val- 
ues uifdc = = 0.263(3) or uogsc — P/v±_ — 0.242(2). However, the good collapses performed 
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with the same 2(7 = 1/2) = 1.45(3) exhibited in the main plots of both figures clearly evidences 
that only longitudinal correlations are relevant in the critical evolution of the system at short times. 

On the other hand, the insets of figured a) and b), show that the finite-size behavior of OP(t) 
is also a power law for the case of large shear field magnitudes, represented by 7 = 10. Opposite to 
the case with 7 = 1/2, the estimated values of tofdc — 0.29(2) and u> gsc = 0.34(3) were in agreement 
with the expected values u>fdc = P/ v \\ = 0.291(2) or ojgsc — P/v±_ = 0.342(2) calculated from 
Table 3. The good collapses performed with z(j = 10) displayed in the main plots of the figures 
allow us to conclude that the same behavior observed for small j's is also exhibited by systems 
with large values of the external fields. 
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Figure 9: Log-log plots of the scaling functions f(t/L z ) (see equation ([HI)) of the time evolution 
of OP{t) when the system is started from FDC's and GSC's configurations (main plot of figure 
a) and b) respectively) corresponding to a system with external field of magnitude 7 = 10. The 
insets of both figures show the finite-size dependence and the corresponding power-law fit in a 
double logarithmic scale. The values of the corresponding fits are ujfdc = 0.29(2) (inset of a) ) and 
Ugsc = 0.34(3) (inset of b)), respectively. 

To summarize, we have shown that the critical dynamic evolution of the system, started from 
either FDC's or GSC initial configurations, can be scaled with the dynamic critical exponent z 
proposed in Section [3J Based on the arguments exposed in Section [3j we can conclude that, in 
the short time limit of the critical evolution, the correlations along the field axis (longitudinal) are 
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more relevant than transverse (perpendicular) correlations. Furthermore, the obtained finite-size 
exponents are only in accordance with those calculated using the critical exponents enlisted in 
Table 3 corresponding to 7 = 10, while at 7 = 1/2 they differ by a factor of nearly 2. 

5 Discussion and Conclusions 

In this work, the second-order phase transition in the 2d non-conserved Ising model under the action 
of an external driving shear field was investigated by studying the critical evolution of the system 
in the short-time regime. In order to apply this method, the dynamic evolution of the system at 
T ~ T c (j) was monitored when it is initiated from fully disordered initial configurations (FDC), 
and from the completely ordered configuration (GSC). 

Starting the system from FDC's configurations, the time evolution of the order parameter OP 
follows a power law behavior at the critical temperature T c , while at slightly different values of T 
the power law is modulated by a scaling function that bends upwards or downwards depending if 
the temperature is less or greater than T c , respectively. The critical evolution was studied on square 
and rectangular lattices of different sizes L x and L y , and the results indicate that the short-time 
critical evolution is free of shape effects. Furthermore, the saturation value reached by OP depends 
only on L x . The critical evolution started from FDC's configurations was studied for different values 
of 7, and the diagram of reduced temperatures T c (7)/T c (0) — 1 versus 7, T c (0) = 2.269 J/ks was 
drawn. As a first observation, all the values found for T c (j) are always greater than the 2d critical 
temperature of the Ising model, that is typical for models driven out of equilibrium by an external 
field. Furthermore, two regimes can be distinguished: 1) a growing regime where T c (j)/T c (0) — 1 
cx 7*. The exponent tp was calculated by means of a linear regression fit, giving -0 = 0.52(3), 
which is consistent with theoretical predictions in ref. [7]; 2) a saturation regime, where T c (j) does 
not change appreciably with 7. In this regime T c (j) ~ 1.18T C (0). A similar diagram was already 
observed in the DLG model [10] , where the temperature grows with the magnitude of the driving 
field and then saturate at large values. 

On the other hand, the critical dynamic behavior of the model was also investigated when it is 
initiated from the ground state configuration (GSC). A decreasing power law is observed for OP(t) 
at the same T c (7) found when the system was started from FDC configurations. This evidences that 
the model experiences a second order phase transition, as expected. Also in this case, the critical 
evolution of the system was simulated on rectangular and square lattices of different sizes. It was 
found that the critical evolution is independent of the shape of the lattice and the saturation value 
of OP only depends on the longitudinal size L x , as for evolutions initiated from FDC configurations. 
So, it is concluded as a general result that in the short-time scale, the system critical evolution is 
free of shape effects, as it is also observed for the critical evolution of the DLG model in the same 
time interval [21]. 

Then, the quantities c, (i = 1,4), defined in sect. 3, were studied in order to calculate the critical 
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exponents of the transition. Starting from FDC's initial configurations, the dynamic critical behav- 
ior at small 7 is slower than at larger 7's. As it can be observed in figure [2] and in Table 1, the order 
parameter exponent c\ is smaller at 7 = 1/32 than at larger 7's. Also, its logarithmic derivative 
exponent C2 is different for this case. On the other hand, the values of c\ and C2 are more stable 
for larger 7's, although c\ at 7 = 50 is slightly smaller than the estimated values for 7 = 5,10. A 
similar scenario was found for the order parameter C3 and its derivative C4 exponents starting from 
the GSC configuration. 

By combining these exponents, the static and dynamic critical exponents /3, is±, v\\ and z were 
calculated and are enlisted in table 3. The order parameter critical exponent for 7 = 1/32 is 
similar to the value calculated for the Ising model, but the values of z, v\\ and v± show that the 
anisotropy introduced by such a small external field is relevant. At large 7's, all the exponents 
are similar within a small range, suggesting that the critical behavior of the model is practically 
independent of the magnitude of the field in this regime. Furthermore, Table 3 also shows that the 
critical exponents of the sheared model do not belong to the universality class of the Ising or of the 
DLG model, even if this model shows a similar phase behavior. 

Finally, the critical exponents summarized above were computed based on the fact that only the lon- 
gitudinal correlation length is relevant for the dynamic critical behavior of the model, independently 
of the initial configuration. In order to check this, a finite-size scaling of the dynamic evolution 
of OP was performed with the aid of equation (fl2|) . This equation must contain both the critical 
dynamic behavior of equations © and (fit)]) at early times of evolution, and also the finite-size 
critical behavior at long times. As a consequence, the finite-size exponents must be ujfdc — @/ v \\ 
and locsc = P/v± for both initial conditions respectively. By measuring the saturated values of 
OP, OP sa t, and computing the exponents ujfdc and uigsc the time series of OP were collapsed 
for the system initiated from ordered and disordered configurations as the main plots of figures [5] 
and [9] show. Therefore, it is concluded that only the longitudinal correlation length £|| takes part 
in the critical evolution of the Ising model when an external shear field is applied. Furthermore, 
the finite-size exponents (jJfdc and uiqsc were not consistent with the rates (3/u« and j3/v± for 
the cases corresponding to 7 = 1/2, while they are in good agreement for a shear field of magni- 
tude 7 = 10. This discrepancy between the predicted and measured critical exponents for the case 
7 = 1/2, together with the fact that the values of the critical exponents estimated for smaller 7 
are not similar with those corresponding to larger values of 7 (see Table [3]) , may be explained by 
conjecturing that, at such small values of the shear rate, the system is less perturbed by the external 
driving. This means that the growth of transverse critical correlations is less affected by the shear 
field, and may become relevant in the short time regime. If this is so, our scaling assumptions will 
be not longer valid, and both £'s need to be considered in order to propose scaling forms for the 
dynamic critical behavior of the model in this regime. In order to study this, new simulations of 
the model with 7 = 1/32 were performed, but they demanded a lot of computational time, specially 
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when the system is started from the GSC configurations because they needed larger lattice sizes 
and evolution time intervals in order to obtain good power laws and saturation regimes (L > 10000, 
evolution time intervals of the order of 10 6 MCS or larger), so we did not obtained reliable results. 
As a consequence, this interesting subject will be left for further research in the future. In contrast, 
at larger j's, the good agreement between /3/v\\ and /3/v±, calculated from the exponents in Table 
[3l and the estimated values of ujfdc and ojgsc respectively, suggest that the critical behavior is 
different from both the cases with smaller j's, and also from the Ising and DLG models at large 
external fields [2T1I22] . 
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